#-------------------------------------------------------------------------------
# Functions to run simulations
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
function gen_ysim(ce::CE_Economy_PInf, T,N; y0=1.)
#-------------------------------------------------------------------------------
    yvals = Array{Float64}(T,N)
    yvals[1,:]  = zeros(N)
    rand_mat    = rand(Normal(0, 1), T,N)
    for nn=1:N
    for tt=2:T
    yvals[tt,nn] = yvals[tt-1,nn]*ce.ρ + ce.η*rand_mat[tt,nn]
    end
    end
    return max.( min.( exp.(yvals), ce.y_grid[end]), ce.y_grid[1])
end

#-------------------------------------------------------------------------------
function gen_type_sim(ce::CE_Economy_PInf, T,N; Initial_Type=1 )
#-------------------------------------------------------------------------------
    type_sim      = zeros(Int, T,N)
    type_sim[1,:] = ones(Int,N).*Initial_Type # Patient
    rand_num      = rand(T,N)
    for nn=1:N
    for tt=2:T
                if      type_sim[tt-1,nn]==1 && rand_num[tt,nn]<ce.ΠT[1,1]  #P-type remains P-type
                        type_sim[tt,nn]=1;
                elseif  type_sim[tt-1,nn]==1 && rand_num[tt,nn]>=ce.ΠT[1,1]  #P-type switches
                        type_sim[tt,nn]=2;
                elseif  type_sim[tt-1,nn]==2 && rand_num[tt,nn]<ce.ΠT[2,2]  #I-type remains I-type
                        type_sim[tt,nn]=2;
                elseif  type_sim[tt-1,nn]==2 && rand_num[tt,nn]>=ce.ΠT[2,2]  #I-type switches
                        type_sim[tt,nn]=1;
                else
                        println("ERROR")
                end
        end
        end
return type_sim
end
#-------------------------------------------------------------------------------



################# ------------------------------------------------------------------------------------- #################
################# ------------------------- *** SIMULATION OF THE ECONOMY *** ------------------------- #################
################# ------------------------------------------------------------------------------------- #################

# Simulation of the economy given realization of all exogenous variables
function simul_econ_PI(ce::CE_Economy_PInf, T,N; T0=2, Initial_Type=1,   Counterfactual=0, Type_Gov="I", b0=0.69, b_pol=0.)
    srand(100)

    #Create TxN exogenous variables
    y_sim         = gen_ysim(ce,  T, N)                                # path for output
    type_sim      = gen_type_sim(ce, T,N, Initial_Type=Initial_Type)   # path for types
    def_rand_sim  = rand(Normal(1.,ce.σ_d),T, N)                       # default shock
    θrand         = rand(T, N)                                         # exit default shocks


#Options for the Argentine Counterfactual - Load Data & Generate Constant Types
#-------------------------------------------------------------------------------
if Counterfactual==1

    # Load Argentine data for the 2006.Q1 - 2012.Q4 period
    cyc_data     =  readdlm(string(Arg_Data_path,Output_data),    ',', Float64)
    y_ARG        =  [1; exp.(cyc_data[:,1]); 1 ];  #first (T0-1) obs are dropped. Last obs is dropped
    T            =  length(y_ARG);
    y_sim        =  repmat(y_ARG,1,N);

    # Set Types for each case [always I-type of always P-type]
    if Type_Gov=="I"
    type_sim     = [ones(5,N); 2*ones(T-5,N)];   #Assume I-type since 2007.Q1
    elseif Type_Gov=="P"
    type_sim     = ones(T,N);     #Assume P-type always
    else
    println(Error)
    end
    def_rand_sim = ones(T,N);

end
#-------------------------------------------------------------------------------


    #Create matrices to fill-in - endogenous variables
    def_sim     = zeros(Int, T, N);
    def_pol     = zeros(Int, T, N);
    b_sim       = zeros(T, N).*NaN;
    π_sim       = zeros(T, N).*NaN;
    q_sim       = zeros(T, N).*NaN;
    SP_sim      = zeros(T, N).*NaN;
    SP_HR_sim   = zeros(T, N).*NaN;
    c_sim       = zeros(T, N).*NaN;
    tb_sim      = zeros(T, N).*NaN;
    def_id      = zeros(T,N);
    VF_sim      = zeros(T, N);

    ## Initialize debt and prior
    b_sim[1:2,:] = b0*ones(2,N);

    ## Interpolations
    knots        = (ce.b_grid, ce.y_grid, 1:2);
    #Prices
    itp_q_mat    = interpolate(knots, ce.q_mat,         Gridded(Linear()));
    #Policies
    itp_def      = interpolate(knots, ce.default,     Gridded(Linear()));
    itp_bp       = interpolate(knots, ce.bp,          Gridded(Linear()));
    #Value Function
    itp_Wf       = interpolate(knots, ce.Wf, Gridded(Linear()));
    itp_Wr       = interpolate(knots, ce.Wr, Gridded(Linear()));
    itp_Wd       = interpolate((ce.y_grid, 1:2), ce.Wd, Gridded(Linear()));


for i=1:N
for t=2:T-1

# Load Current State
Ti, b_v, y_v = type_sim[t,i], b_sim[t,i], y_sim[t,i];

VF_sim[t,i]  = itp_Wf[b_v, y_v, Ti]


# IF COUNTRY IS NOT IN A DEFAULT OR IF IT EXITS THE DEFAULT TODAY
if def_sim[t-1,i]==0 || (def_sim[t-1,i]==1 && θrand[t,i] <= ce.θ)

        # CURRENT DEFAULT CHOICE
        def_v     = (itp_Wr[b_v,y_v,Ti]./itp_Wd[y_v,Ti])

        if def_v.<=def_rand_sim[t,i] # No default today

                    def_sim[t,i] = 0.;
                    def_id[t,i]  = def_id[t-1,i]

                    π_sim[t,i]   = ifelse(Ti==1, 0.0, ce.π_min)

                    ## Bond Policies [same for the 2 types]
                    bp = itp_bp[b_v, y_v, Ti];
                    b_sim[t+1,i] = bp;

                    #Bond policy under the P-type counterfactual
                    if (Counterfactual==1  && size(b_pol,1)>2)
                        bp = b_pol[t,i]
                        b_sim[t+1,i] = bp
                    end

                    # Prices (depend on bp, y, ζ_prime)
                    q_sim[t,i]      = itp_q_mat[b_sim[t+1,i], y_sim[t,i], Ti]
                    q_HR            = itp_q_mat[b_sim[t+1,i], y_sim[t,i], 1]
                    SP_sim[t,i]     = Spreads_fx(ce,q_sim[t,i])    # annualized, pp
                    SP_HR_sim[t,i]  = Spreads_fx(ce,q_HR)          # annualized, pp

                    # Consumption and trade balance
                    c_sim[t,i]   = y_sim[t,i] - b_sim[t,i]* ((1-ce.mb)*ce.zb + ce.mb) - ce.Bval*(1+π_sim[t,i]) + q_sim[t,i]*(b_sim[t+1,i]-(1-ce.mb)*b_sim[t,i])
                    tb_sim[t,i]  = (y_sim[t,i] - c_sim[t,i])     # Y = C + NX


        elseif def_v.>def_rand_sim[t,i] # defaults today

                    def_id[t,i]  = 1 + def_id[t-1,i]
                    def_sim[t,i] = 1;
                    def_pol[t,i] = 1;   #Only relevant to compute default frequency

                    #Output Cost
                    yIdf = y_sim[t,i]-max((ce.d0+ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)
                    yPdf = y_sim[t,i]-max((ce.d0-ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)

                    y_sim[t,i]    = yPdf*(Ti==1) + yIdf*(Ti==2)
                    c_sim[t,i]    = y_sim[t,i]
                    b_sim[t+1,i]  = 0.0
                    π_sim[t,i]    = NaN
                    tb_sim[t,i]   = NaN

          else
          println("ERROR")
          end

# IF COUNTRY IS ALREADY IN A DEFAULT AND IT DOES NOT EXIT TODAY
elseif (def_sim[t-1,i]==1 && θrand[t,i] > ce.θ)

                    def_sim[t,i] = 1
                    def_id[t,i]  = def_id[t-1,i]

                    yIdf = y_sim[t,i]-max((ce.d0+ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)
                    yPdf = y_sim[t,i]-max((ce.d0-ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)

                    y_sim[t,i]    = yPdf*(Ti==1) + yIdf*(Ti==2)
                    c_sim[t,i]    = y_sim[t,i]
                    b_sim[t+1,i]  = 0.0
                    π_sim[t,i]    = NaN
                    tb_sim[t,i]   = NaN


else
println("ERROR")
end

end # for T
end # for N

return  y_sim[T0:end-1,:], type_sim[T0:end-1,:], b_sim[T0:end-1,:], π_sim[T0:end-1,:], def_sim[T0:end-1,:], def_pol[T0:end-1,:], def_id[T0:end-1,:], q_sim[T0:end-1,:],
        SP_sim[T0:end-1,:], SP_HR_sim[T0:end-1,:], c_sim[T0:end-1,:], tb_sim[T0:end-1,:], b_sim[:,:],  VF_sim[T0:end-1,:]
end
